$ontext
Loop version with multiple scenarios for 4 variables. n1=8, n2=1, n3=2, n4=1.

This is DICE-GIS model. We use includes to keep modeling clearn.
It is prepared for published version of model in PNAS.
It needs the four include files to run.

Note that version has the loops. It contains n1 x n2 x n3 x n4 runs and will
be very slow if all scenarios are calculated.

Prepared on March 29, 2019
$offtext

$title        DICE-GIS-March 2019:DICE-GIS-v32-standard-1111-loop.gms
set        t  Time periods (5 years per period)                    /1*500/
** Random sets
set   kkdisc alternative discount rates /1*8/
set   kkdamgis alternative damages GIS /1*1/
set   kkmelt alternative melt rates /1*2/
set   kkvol alternative volume constraints /1*1/
PARAMETERS
** If optimal control
        ifopt    Indicator where optimized is 1 and base is 0      /1/
** Damages if on
       ifdamgis  If damages from GIS melt                         /1/
       ifstandam     If normal damages on                         /1/
** coef on initial gis volume
       coefvolinit  Coefficient initial volume                    /1/
** Availability of fossil fuels
      fosslim  Maximum cumulative extraction fossil fuels (GtC)  /6000/
** Preferences
*     elasmu   Elasticity of marginal utility of consumption     /1.45 /
     elasmu   Elasticity of marginal utility of consumption     /.01 /
     prstp    Initial rate of social time preference per year   /.015  /
*      elasmu   Elasticity of marginal utility of consumption     /.001 /
*      prstp    Initial rate of social time preference per year   /.01  /

** Population and technology
        gama     Capital elasticity in production function        /.300    /
        pop0     Initial world population 2015 (millions)         /7403    /
        popadj   Growth rate to calibrate to 2050 pop projection  /0.134   /
        popasym  Asymptotic population (millions)                 /11500   /
        dk       Depreciation rate on capital (per year)          /.100    /
        q0       Initial world gross output 2015 (trill 2010 USD) /105.5   /
        k0       Initial capital value 2015 (trill 2010 USD)      /223     /
        a0       Initial level of total factor productivity       /5.115    /
        ga0      Initial growth rate for TFP per 5 years          /0.076   /
        dela     Decline rate of TFP per 5 years                  /0.005   /
** Emissions parameters
        gsigma1  Initial growth of sigma (per year)                   /-0.0152 /
        dsig     Decline rate of decarbonization (per period)         /-0.001  /
        eland0   Carbon emissions from land 2015 (GtCO2 per year)     / 2.6    /
        deland   Decline rate of land emissions (per period)          / .115   /
        e0       Industrial emissions 2015 (GtCO2 per year)           /35.85    /
        miu0     Initial emissions control rate for base case 2015    /.03     /

** Carbon cycle
* Initial Conditions
        mat0   Initial Concentration in atmosphere 2015 (GtC)        /851    /
        mu0    Initial Concentration in upper strata 2015 (GtC)      /460    /
        ml0    Initial Concentration in lower strata 2015 (GtC)      /1740   /
        mateq  Equilibrium concentration atmosphere  (GtC)           /588    /
        mueq   Equilibrium concentration in upper strata (GtC)       /360    /
        mleq   Equilibrium concentration in lower strata (GtC)       /1720   /
* Flow paramaters
        b12      Carbon cycle transition matrix                      /.12   /
        b23      Carbon cycle transition matrix                      /0.007 /
* These are for declaration and are defined later
        b11      Carbon cycle transition matrix
        b21      Carbon cycle transition matrix
        b22      Carbon cycle transition matrix
        b32      Carbon cycle transition matrix
        b33      Carbon cycle transition matrix
        sig0     Carbon intensity 2010 (kgCO2 per output 2005 USD 2010)
** Climate model parameters
        t2xco2   Equilibrium temp impact (oC per doubling CO2)    / 3.1  /
        fex0     2015 forcings of non-CO2 GHG (Wm-2)              / 0.5  /
        fex1     2100 forcings of non-CO2 GHG (Wm-2)              / 1.0  /
        tocean0  Initial lower stratum temp change (C from 1900)  /.0068 /
        tatm0    Initial atmospheric temp change (C from 1900)    /0.85  /
        c1       Climate equation coefficient for upper level     /0.1005  /
        c3       Transfer coefficient upper to lower stratum      /0.088   /
        c4       Transfer coefficient for lower level             /0.025   /
        fco22x   Forcings of equilibrium CO2 doubling (Wm-2)      /3.6813  /
** Climate damage parameters
        a10       Initial damage intercept                         /0       /
        a20       Initial damage quadratic term
        a1        Damage intercept                                 /0       /
       a2        Damage quadratic term                            /0.00236 /
        a3        Damage exponent                                  /2.00    /
** Abatement cost
        expcost2  Exponent of control cost function               / 2.6  /
        pback     Cost of backstop 2010$ per tCO2 2015            / 550  /
        gback     Initial cost decline backstop cost per period   / .025 /
        tnopol    Period before which no emissions controls base  / 60   /
        tzeroemis Period of zero emissions                        / 500   /
        cprice0   Initial base carbon price (2010$ per tCO2)      / 2    /
        gcprice   Growth rate of base carbon price per year       /.00002   /
** Long run emissions control rates
        tnonegem1 Time period for no negative emissions 1          / 20    /
        tnonegem2 Time period for no negative emissions 2          / 30    /
        limmiu1    Upper limit on control rate after TNONEGEM1      / 1.2 /
        limmiu2    Upper limit on control rate after TNONEGEM2      / 1. /

** Scaling and inessential parameters.They ensure that MU of first period's consumption =1 and PV cons = PV utilty
        scale1      Multiplicative scaling coefficient           /0.0302455265681763 /
        scale2      Additive scaling coefficient                 /-10993.704/ ;

PARAMETERS
* Loop vectors
        loopdamcoefgis(kkdamgis)
        loopdisc(kkdisc)
        loopelast(kkdisc)
        loopmeltcoef(kkmelt)
        loopvol(kkvol)
        loopscc(kkdamgis,kkmelt,kkdisc,kkvol,t)
        looptemp(kkdamgis,kkmelt,kkdisc,kkvol,t)
        loopmat(kkdamgis,kkmelt,kkdisc,kkvol,t)
        loopdamf(kkdamgis,kkmelt,kkdisc,kkvol,t)
        loopcca(kkdamgis,kkmelt,kkdisc,kkvol,t)
        loopy(kkdamgis,kkmelt,kkdisc,kkvol,t)
        loopem(kkdamgis,kkmelt,kkdisc,kkvol,t)
        loopr(kkdamgis,kkmelt,kkdisc,kkvol,t)
        loopmiu(kkdamgis,kkmelt,kkdisc,kkvol,t)
        loopvolgis(kkdamgis,kkmelt,kkdisc,kkvol,t)
        loopobjective(kkdamgis,kkmelt,kkdisc,kkvol)
        loopvdot(kkdamgis,kkmelt,kkdisc,kkvol,t)
        looptstar(kkdamgis,kkmelt,kkdisc,kkvol,t)
        loopsign(kkdamgis,kkmelt,kkdisc,kkvol,t)
        looptd(kkdamgis,kkmelt,kkdisc,kkvol,t)
        scc(t)
;

parameter  loopdamcoefgis(kkdamgis)  Coefficients GIS damage
/  1  0.0002   / ;
parameter loopmeltcoef(kkmelt) Coefficients for melt rate multiplier
/  1  1, 2   2  /  ;
parameter  loopvol(kkvol)  Coefficients  for  minimum volume
/  1  1  /  ;
 parameter loopdisc(kkdisc) Coefficients for pure rate time pref
/  1  .015, 2 .05, 3 .04, 4 .03, 5 .02, 6 .01, 7 .001 , 8 .001/
 parameter loopelast(kkdisc) Coefficients for pure rate time pref
/  1  1.45 ,2 .01 , 3 .01, 4 .01, 5 .01, 6 .01, 7 .99, 8 .01 /

;
* Program control variables
sets     tfirst(t), tlast(t), tearly(t), tlate(t);

$include GIS_inc_1_v32

PARAMETERS
        l(t)          Level of population and labor
        al(t)         Level of total factor productivity
        sigma(t)      CO2-equivalent-emissions output ratio
        rr(t)         Average utility social discount rate
        ga(t)         Growth rate of productivity from
        forcoth(t)    Exogenous forcing for other greenhouse gases
        gl(t)         Growth rate of labor
        gcost1        Growth of cost factor
        gsig(t)       Change in sigma (cumulative improvement of energy efficiency)
        etree(t)      Emissions from deforestation
        cumetree(t)   Cumulative from land
        cost1(t)      Adjusted cost for backstop
        lam           Climate model parameter
        gfacpop(t)    Growth factor population
        pbacktime(t)  Backstop price
        optlrsav      Optimal long-run savings rate used for transversality
        scc(t)        Social cost of carbon
        photel(t)     Carbon Price under no damages (Hotelling rent condition)
        ppm(t)        Atmospheric concentrations parts per million
        atfrac(t)     Atmospheric share since 1850
        atfrac2010(t)     Atmospheric share since 2010 ;
* Program control definitions
        tfirst(t) = yes$(t.val eq 1);
        tlast(t)  = yes$(t.val eq card(t));
* Parameters for long-run consistency of carbon cycle
        b11 = 1 - b12;
        b21 = b12*MATEQ/MUEQ;
        b22 = 1 - b21 - b23;
        b32 = b23*mueq/mleq;
        b33 = 1 - b32 ;
* Further definitions of parameters
        a20 = a2;
        sig0 = e0/(q0*(1-miu0));
        lam = fco22x/ t2xco2;
        l("1") = pop0;
        loop(t, l(t+1)=l(t););
        loop(t, l(t+1)=l(t)*(popasym/L(t))**popadj ;);

        ga(t)=ga0*exp(-dela*5*((t.val-1)));
        al("1") = a0; loop(t, al(t+1)=al(t)/((1-ga(t))););
        gsig("1")=gsigma1; loop(t,gsig(t+1)=gsig(t)*((1+dsig)**5) ;);
        sigma("1")=sig0;   loop(t,sigma(t+1)=(sigma(t)*exp(gsig(t)*5)););

        pbacktime(t)=pback*(1-gback)**(t.val-1);
        cost1(t) = pbacktime(t)*sigma(t)/expcost2/1000;
        etree(t) = eland0*(1-deland)**(t.val-1);
        cumetree("1")= 100; loop(t,cumetree(t+1)=cumetree(t)+etree(t)*(5/3.666););
        rr(t) = 1/((1+prstp)**(5*(t.val-1)));
        forcoth(t) = fex0+ (1/17)*(fex1-fex0)*(t.val-1)$(t.val lt 18)+ (fex1-fex0)$(t.val ge 18);
        optlrsav = (dk + .004)/(dk + .004*elasmu + prstp)*gama;

*Freeze sigma for high T for neg emissions
 sigma(t)$(t.val gt 40) = sigma("40");

VARIABLES
        MIU(t)          Emission control rate GHGs
        FORC(t)         Increase in radiative forcing (watts per m2 from 1900)
        TATM(t)         Increase temperature of atmosphere (degrees C from 1900)
        TOCEAN(t)       Increase temperatureof lower oceans (degrees C from 1900)
        MAT(t)          Carbon concentration increase in atmosphere (GtC from 1750)
        MU(t)           Carbon concentration increase in shallow oceans (GtC from 1750)
        ML(t)           Carbon concentration increase in lower oceans (GtC from 1750)
        E(t)            Total CO2 emissions (GtCO2 per year)
        EIND(t)         Industrial emissions (GtCO2 per year)
        C(t)            Consumption (trillions 2005 US dollars per year)
        K(t)            Capital stock (trillions 2005 US dollars)
        CPC(t)          Per capita consumption (thousands 2005 USD per year)
        I(t)            Investment (trillions 2005 USD per year)
        S(t)            Gross savings rate as fraction of gross world product
        RI(t)           Real interest rate (per annum)
        Y(t)            Gross world product net of abatement and damages (trillions 2005 USD per year)
        YGROSS(t)       Gross world product GROSS of abatement and damages (trillions 2005 USD per year)
        YNET(t)         Output net of damages equation (trillions 2005 USD per year)
        DAMAGES(t)      Damages (trillions 2005 USD per year)
        DAMFRAC(t)      Damages as fraction of gross output
        ABATECOST(t)    Cost of emissions reductions  (trillions 2005 USD per year)
        MCABATE(t)      Marginal cost of abatement (2005$ per ton CO2)
        CCA(t)          Cumulative industrial carbon emissions (GTC)
        CCATOT(t)       Total carbon emissions (GtC)
        PERIODU(t)      One period utility function
        CPRICE(t)       Carbon price (2005$ per ton of CO2)
        CEMUTOTPER(t)   Period utility
        UTILITY         Welfare function;

NONNEGATIVE VARIABLES  MIU, MAT, MU, ML, Y, YGROSS, C, K, I;

EQUATIONS
*Emissions and Damages
        EEQ(t)           Emissions equation
        EINDEQ(t)        Industrial emissions
        CCACCA(t)        Cumulative industrial carbon emissions
        CCATOTEQ(t)        Cumulative total carbon emissions
        FORCE(t)         Radiative forcing equation
        DAMFRACEQ(t)     Equation for damage fraction
        DAMEQ(t)         Damage equation
        ABATEEQ(t)       Cost of emissions reductions equation
        MCABATEEQ(t)     Equation for MC abatement
        CARBPRICEEQ(t)   Carbon price equation from abatement

*Climate and carbon cycle
        MMAT(t)          Atmospheric concentration equation
        MMU(t)           Shallow ocean concentration
        MML(t)           Lower ocean concentration
        TATMEQ(t)        Temperature-climate equation for atmosphere
        TOCEANEQ(t)      Temperature-climate equation for lower oceans

*Economic variables
        YGROSSEQ(t)      Output gross equation
        YNETEQ(t)        Output net of damages equation
        YY(t)            Output net equation
        CC(t)            Consumption equation
        CPCE(t)          Per capita consumption definition
        SEQ(t)           Savings rate equation
        KK(t)            Capital balance equation
        RIEQ(t)          Interest rate equation

* Utility
        CEMUTOTPEREQ(t)  Period utility
        PERIODUEQ(t)     Instantaneous utility function equation
        UTIL             Objective function      ;

$include GIS_inc_2_v32

** Equations of the model
*Emissions and Damages
 eeq(t)..             E(t)           =E= EIND(t) + etree(t);
 eindeq(t)..          EIND(t)        =E= sigma(t) * YGROSS(t) * (1-(MIU(t)));
 ccacca(t+1)..        CCA(t+1)       =E= CCA(t)+ EIND(t)*5/3.666;
 ccatoteq(t)..        CCATOT(t)      =E= CCA(t)+cumetree(t);
 force(t)..           FORC(t)        =E= fco22x * ((log((MAT(t)/588.000))/log(2))) + forcoth(t);
 dameq(t)..           DAMAGES(t)     =E= YGROSS(t) * DAMFRAC(t);
 abateeq(t)..         ABATECOST(t)   =E= YGROSS(t) * cost1(t) * (MIU(t)**expcost2);
 mcabateeq(t)..       MCABATE(t)     =E= pbacktime(t) * MIU(t)**(expcost2-1);
 carbpriceeq(t)..     CPRICE(t)      =E= pbacktime(t) * (MIU(t))**(expcost2-1);

*Substitute equation from standard DICE
 damfraceq(t) ..      DAMFRAC(t)     =E= ifstandam*(a1*TATM(t)+a2*TATM(t)*TATM(t))+ifdamgis*cgisdam*(volzero-GISVOL(T)) ;

*Climate and carbon cycle
 mmat(t+1)..          MAT(t+1)       =E= MAT(t)*b11 + MU(t)*b21 + (E(t)*(5/3.666));
 mml(t+1)..           ML(t+1)        =E= ML(t)*b33  + MU(t)*b23;
 mmu(t+1)..           MU(t+1)        =E= MAT(t)*b12 + MU(t)*b22 + ML(t)*b32;
 tatmeq(t+1)..        TATM(t+1)      =E=  TATM(t) + c1 * ((FORC(t+1)-(fco22x/t2xco2)*TATM(t))-(c3*(TATM(t)-TOCEAN(t))));

* tatmeq(t+1)..        TATM(t+1)      =E= 6$(t.val lt 100)+0.$(t.val ge 100);
toceaneq(t+1)..      TOCEAN(t+1)    =E= TOCEAN(t) + c4*(TATM(t)-TOCEAN(t));

*Economic variables
 ygrosseq(t)..        YGROSS(t)      =E= (al(t)*(L(t)/1000)**(1-GAMA))*(K(t)**GAMA);
 yneteq(t)..          YNET(t)        =E= YGROSS(t)*(1-damfrac(t));
 yy(t)..              Y(t)           =E= YNET(t) - ABATECOST(t);
 cc(t)..              C(t)           =E= Y(t) - I(t);
 cpce(t)..            CPC(t)         =E= 1000 * C(t) / L(t);
 seq(t)..             I(t)           =E= S(t) * Y(t);
 kk(t+1)..            K(t+1)         =E= (1-dk)**5 * K(t) + 5 * I(t);
 rieq(t+1)..          RI(t)          =E= (1+prstp) * (CPC(t+1)/CPC(t))**(elasmu/5) - 1;

*Utility
 cemutotpereq(t)..    CEMUTOTPER(t)  =E= PERIODU(t) * L(t) * rr(t);
 periodueq(t)..       PERIODU(t)     =E= ((C(T)*1000/L(T))**(1-elasmu)-1)/(1-elasmu)-1;
 util..               UTILITY        =E= -gisvol('2')*coefvolinit+5 * scale1 * sum(t,  CEMUTOTPER(t)) + scale2 ;

*Resource limit
 CCA.up(t)       = fosslim;

* Control rate limits

 MIU.up(t) = limmiu1;
 MIU.up(t)$(t.val<tnonegem1) = 1;
 MIU.up(t)$(t.val>tnonegem2) = limmiu2;

**  Upper and lower bounds for stability
K.LO(t)         = 1;
MAT.LO(t)       = 10;
MU.LO(t)        = 100;
ML.LO(t)        = 1000;
C.LO(t)         = 2;
TOCEAN.UP(t)    = 20;
TOCEAN.LO(t)    = -1;
TATM.UP(t)      = 20;
CPC.LO(t)       = .01;

* Control variables
s.fx(t)$(t.val gt 100) = .25;

* Initial conditions
CCA.FX(tfirst)    = 400;
K.FX(tfirst)      = k0;
MAT.FX(tfirst)    = mat0;
MU.FX(tfirst)     = mu0;
ML.FX(tfirst)     = ml0;
TATM.FX(tfirst)   = tatm0;
TOCEAN.FX(tfirst) = tocean0;

$include GIS_inc_3_v32

MIU.lo(t)$(t.val>100) = 1.0;
** Solution options
option iterlim = 99900;
option reslim = 99999;
option solprint = on;
option limrow = 0;
option limcol = 0;
model  CO2 /all/;

loop(kkmelt, loop(kkvol, loop(kkdisc, loop(kkdamgis,

** Alternative parameters
meltmult=loopmeltcoef(kkmelt);
avoldot=avoldot0*meltmult;
prstp= loopdisc(kkdisc);
elasmu= loopelast(kkdisc);
rr(t) = 1/((1+prstp)**(5*(t.val-1)));
cgisdam = loopdamcoefgis(kkdamgis);
gisvol.lo(t)=loopvol(kkvol);

GISVOL.fx(tfirst) = volzero;

* For base run, this subroutine calculates Hotelling rents
If (ifopt eq 0,
       a2 = 0.0001;
       cgisdam = .000001;
       solve CO2 maximizing UTILITY using nlp;
       solve CO2 maximizing UTILITY using nlp;
       solve CO2 maximizing UTILITY using nlp;
       photel(t)=cprice.l(t);
       a2 = a20;
       cgisdam = loopdamcoefgis(kkdamgis);
      cprice.up(t)$(t.val<tnopol) = photel(t);
);

solve co2 maximizing utility using nlp;
solve co2 maximizing utility using nlp;
solve co2 maximizing utility using nlp;

** POST-SOLVE. Calculate social cost of carbon and other variables
scc(t) = -1000*eeq.m(t)/(.00001+cc.m(t));
atfrac(t) = ((mat.l(t)-588)/(ccatot.l(t)+.000001  ));
atfrac2010(t) = ((mat.l(t)-mat0)/(.00001+ccatot.l(t)-ccatot.l('1')  ));
ppm(t)    = mat.l(t)/2.13;
* Loop variables
loopscc(kkdamgis,kkmelt,kkdisc,kkvol,t) = -1000*eeq.m(t)/(.00000001+cc.m(t));
looptemp(kkdamgis,kkmelt,kkdisc,kkvol,t) = tatm.l(t);
loopmat(kkdamgis,kkmelt,kkdisc,kkvol,t) = mat.l(t);
loopy(kkdamgis,kkmelt,kkdisc,kkvol,t) = y.l(t);
loopem(kkdamgis,kkmelt,kkdisc,kkvol,t) = eind.l(t);
loopr(kkdamgis,kkmelt,kkdisc,kkvol,t) = ri.l(t);
loopmiu(kkdamgis,kkmelt,kkdisc,kkvol,t) = miu.l(t);
loopvolgis(kkdamgis,kkmelt,kkdisc,kkvol,t) = gisvol.l(t);
loopdamf(kkdamgis,kkmelt,kkdisc,kkvol,t) = damfrac.l(t);
loopcca(kkdamgis,kkmelt,kkdisc,kkvol,t) = cca.l(t);
loopobjective(kkdamgis,kkmelt,kkdisc,kkvol) = utility.l;
 loopvdot(kkdamgis,kkmelt,kkdisc,kkvol,t)  = vdot.l(t);
 looptstar(kkdamgis,kkmelt,kkdisc,kkvol,t)  = tstarlin.l(t);
 looptd(kkdamgis,kkmelt,kkdisc,kkvol,t)  = td.l(t);
scc(t) = -1000*eeq.m(t)/(.00001+cc.m(t));
 loopsign(kkdamgis,kkmelt,kkdisc,kkvol,t)  = signtd.l(t);
););););

* Produces an output file.
* The statement at the end of the *.lst file "Output..." will tell you where to find the file.
file results /DICE-GIS-v32-standard-1111-loop.csv/;     results.nd = 10 ; results.nw = 0 ; results.pw=20000; results.pc=5;
put results;
put /"Results of DICE-GIS-v32-standard-1111-loop.gms";
$include GIS_inc4_v32

display s.l, cprice.l, miu.l, cprice.l,scc,ri.l, gisvol.l;
display Loopdamcoefgis,  loopmeltcoef,  loopvol,  loopdisc, loopelast;
display scc,tstarlin.l,td.l,signtd.l,tatm.l;


